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Abstract 

The study makes use of polynomial chaos expansions to compute 
Sobol’ indices within the frame of a global sensitivity analysis of hydro- 
dispersive parameters in a simplified vertical cross-section of a segment 
of the subsurface of the Paris Basin. Applying conservative ranges, 
the uncertainty in 78 input variables is propagated upon the mean life¬ 
time expectancy of water molecules departing from a specific location 
within a highly confining layer situated in the middle of the model do¬ 
main. Lifetime expectancy is a hydrogeological performance measure 
pertinent to safety analysis with respect to subsurface contaminants, 
such as radionuclides. The sensitivity analysis indicates that the vari¬ 
ability in the mean lifetime expectancy can be sufficiently explained 
by the uncertainty in the petrofacies, i.e. the sets of porosity and hy¬ 
draulic conductivity, of only a few layers of the model. The obtained 
results provide guidance regarding the uncertainty modeling in future 
investigations employing detailed numerical models of the subsurface 
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of the Paris Basin. Moreover, the study demonstrates the high ef¬ 
ficiency of sparse polynomial chaos expansions in computing Sobol’ 
indices for high-dimensional models. 

Keywords - global sensitivity analysis - polynomial chaos ex¬ 
pansions - groundwater flow - lifetime expectancy - deep geological 
storage 


1 Introduction 


With the improvement of computing power, numerical modeling has become 
a popular tool for understanding and predicting various kinds of subsurface 
processes addressed in the fields of geology and hydrogeology. However, the 
incomplete/imprecise knowledge of the underground system frequently com¬ 
pels the modeller to make a number of approximations and assumptions with 
regard to the geometry of geological structures, the presence of discontinu¬ 
ities and/or the spatial distribution of hydro-dispersive parameters in their 


models Renard (2007). These uncertainties can possibly lead to large vari¬ 


abilities in the predictive modeling of subsurface processes and thus, it be¬ 
comes of major importance to account for the aforementioned assumptions 
in the frame of uncertainty and sensitivity analyses. Uncertainty analysis 
(UA) aims at quantifying the variability of a given response of interest as a 
function of uncertain input factors, whereas sensitivity analysis (SA) has the 
purpose to identify the input factors responsible for this variability. Hence, 
SA determines the key variables to be described in further detail in order to 
reduce the uncertainty in the predictions of a model. 

Methods of SA are typically classified in two categories: local SA and 
global SA methods. The former investigate effects of variations of the input 
factors in the vicinity of nominal values, whereas the latter aim at quanti¬ 
fying the output uncertainty due to variations in the input factors in their 
entire domain. Among several global SA methods proposed in the litera¬ 
ture, of interest herein is SA with Sobol’ sensitivity indices , which belongs 
to the broader class of variance-based methods Saltelli et al. (2008). These 


methods rely upon the decomposition of the response variance as a sum of 
contributions of each input factor or combinations thereof and do not assume 
any kind of linearity or monotonicity of the model. We note that the Fourier 


(1999) 


amplitude sensitivity test (FAST) Cukier et al. (1978); Saltelli et al. 
indices enter this class as well. 

Various methods have been investigated for computing the Sobol’ indices 
that were first defined in Sobol’ (1993), see e.g. Archer et al. (1997); Sobol’ 
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(2001); Saltelli (2002); Sobol’ and Kucherenko (2005); Saltelli et al. (2010). 


In these papers, Monte Carlo simulation is used as a tool to estimate these 
sensitivity indices. This has revealed extremely costly, although more ef¬ 


ficient estimators have been recently proposed Sobol et al. (2007); Janon 


et al. (2013). In the recent years, new approaches using surrogate models 


have been introduced in the field of global SA Oakley and O’Hagan (2004); 


Marrel et al. (2009); Storlie et al. (2009); Zuniga et al. (2013). A popu¬ 


lar method to compute the Sobol’ indices, originally introduced by |Sudret 


(2008), is by post-processing the coefficients of a polynomial chaos expansion 
(PCE) meta-model of the response quantity of interest. PCE constitutes an 
efficient UA method in which the key concept is to expand the model re¬ 
sponse onto a basis made of orthogonal polynomials in the input variables. 
Once a PCE representation is available, the Sobol’ indices can be calculated 
analytically with elementary operations at almost no additional computa¬ 
tional cost. Sparse PCE make the approach even more efficient, as shown in 


Blatman and Sudret (2010a). 


In the frame of the stochastic modeling of subsurface flow and mass trans¬ 
port, PCE meta-models have proven to be comprehensive and robust tools for 
performing SA at low computational cost. As an example, applying a PCE- 
based global SA upon a fine-grid numerical model of flow and mass transport 


in a heterogeneous porous medium, Fajraoui et al. (2011) and Younes et al. 


(2013) established the transient effect of uncertain flow boundary conditions, 
hydraulic conductivities and dispersivities on solute concentrations at given 
observation points. Sochala and Le Maitre (2013) propagated uncertain soil 
parameters upon three different physical models of subsurface unsaturated 
flow. Their study proved the higher efficiency of PCE meta-models, in com¬ 
parison to a classical Monte-Carlo method, for representing the variability of 
the output quantity at low computational cost. In the frame of radionuclide 
transport simulation in aquifers, Ciricllo et al. (2013) analyzed the statistical 
moments of the peak solute concentration measured at a specific location, as 
a function of the conductivity field, the dispersivity coefficients and the par¬ 
tition coefficients associated to the heterogeneous media. The comparison of 
the Sobol’ indices obtained for various degrees of PCE meta-models showed 
that low-degree models can yield reliable indices while considerably reducing 
the computational burden. Fornraggia et al. (2013) used PCE-based sensi¬ 


tivity indices to investigate effects of uncertainty in hydrogeological variables 
on the evolution of a basin-scale sedimentation process. However, the various 
aforementioned contributions consider simplified models for the description 
of subsurface flow and mass transport. A detailed site characterization model 
was employed by Laloy et al. (2013), but global SA was confined to flow pro¬ 
cesses. 
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In the scope of the deep geological storage of radioactive wastes, ANDRA 
(French National Radioactive Waste Management Agency) has conducted 
several studies to assess the potentiality of a clay-rich layer for establish¬ 
ing a mid to long-lived radioactive waste disposal in the subsurface of the 
Paris Basin. The thick impermeable layer from Callovo-Oxfordian (COX) age 


has been extensively studied (Delay et ah, 2006; Distinguin and Lavanchy 


2007; Enssle et ah, 2011) together with the two major limestone aquifers, 


in place of the Dogger and the Oxfordian sequences (Brigaud et al. 2010 


Linard et al. 

2011 

Landrein et al., 

2013) 

encompassing the claystone forma- 

tion. A recent study ( 

Deman et al. 

20L 

5) used a high-resolution integrated 


Meuse/Haute-Marne hydrogeological model (AND, 2012) to compute the av¬ 
erage time for water molecules departing from a given area in the COX to 
reach the limits of the domain where the numerical model is defined. SA 
over hydro-dispersive parameters in 14 hydrogeological layers proved that 
the Dogger and Oxfordian limestone sequences have a large influence on the 
residence time of groundwater. Indeed, advection processes occurring in per¬ 
meable layers strongly influence the water transit in the subsurface of the 
Paris Basin, in contrast to the slow-motion diffusive processes taking place 
in impermeable rocks. 

However, the analysis of the effect of uncertainties related to other advective- 
dispersive parameters, such as boundary conditions, orientations and anisotropies 
of hydraulic conductivity tensors or magnitudes of dispersion parameters, 
represents a great effort that cannot be carried out with the integrated model 
at reasonable computational costs. Addressing the issue of performing UA 
with the use of high-resolution numerical models of geological reservoirs, 
Castellini and co-workers (Castellini et al. 2003) established that numerical 
models built at the coarse scale, but covering a reasonable number of geolog¬ 
ical and geostatistical features, can be particularly informative in capturing 
the main subsurface processes at low computational costs. 

The present study introduces a vertical two-dimensional multi-layered hy¬ 
drogeological model representing a simplification of the underground media 
of the Paris Basin in the vicinity of the site of Bure and does not integrate the 
complex geometry of the layers, neither does it include the numerous discon¬ 
tinuities or heterogeneities observed in the held. It must be emphasized that 
this simplified model is not aimed at site characterization, but at evaluation 
of hydrogeological performance through SA for calibration purposes. 

The main objective of the present work is to assess the effect of multiple 
advective-dispersive parameters on the mean lifetime expectancy (MLE) of 
water molecules departing from a target zone in the central layer. The MLE 
corresponds to the average time required for a given solute at a specific lo¬ 
cation to reach any outlet of the model domain and is a critical quantity in 
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safety analysis dealing with subsurface contaminants such as radionuclides. 


This work represents a substantial complement to the study by Deman etah 


(2015) by encompassing a large scope of uncertain factors, which cannot 


be assessed using the integrated model due to the computational burden. 
Conservative uncertainty ranges are defined for the input factors analyzed 
in the frame of a SA relying on the estimation of PCE-based Sobol’ in¬ 
dices. The sparse PCE approach was chosen because of its ability to tackle 
high-dimensional problems with great efficiency. The study provides recom¬ 
mendations for future investigations employing the high-resolution integrated 
Meuse/Haute-Marne hydrogeological numerical model of the Paris Basin; in 
particular, it identifies the sets of parameters that can be fixed to their nom¬ 
inal values without significantly affecting the MLE variability as well as the 
sets of parameters to be described in further detail. 

The paper is organized as follows: The subsequent Section 2 provides 
a comprehensive description of the considered hydrogeological model, while 
Section 3 presents the concepts of SA with Sobol’ indices and the compu¬ 
tation of those indices using PCE. Section 4 includes the results of UA and 
SA of the model, along with interpretations accounting for the underlying 
physics. Finally, the conclusive Section 5 summarizes the study, highlighting 
the main findings of the above analyses, and provides recommendations for 
future investigations. 

2 The numerical model 


2.1 Geometry and finite element mesh 


Originally inspired by the COUPLEX numerical model from Bourgeat et ah 


(2004), the present model stands as a vertical two-dimensional (x-z) cross- 
section of 25,000 x 1,040 meters representing a segment of the Paris Basin 
subsurface. The mesh is discretized into 5x5 meters square elements for 
a total of 1,040,000 elements. In order to subdivide the domain into en¬ 
tities related to geological formations, the main features of the subsurface 
were extracted from the lithostratigraphic log of the deep EST433 borehole 


(Landrein et ah, 2013) in the vicinity of the experimental site of Bure (Haute- 


Marne, France). Therefore, the model consists of 15 hydrogeological layers 
characterized by tabular geometries, uniform thicknesses and homogeneous 
parameters. Figure [T| summarizes the geometry of the model and gives an 
overview on the succession and thicknesses of layers. 

The bottom layer stands as a 110 m thick low-permeability layer at¬ 
tributed to the Toarcian marl formation (T). Overlying the latter, the suc- 
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cession of carbonate formations from the Dogger sequence is subdivided into 
5 layers of which the total thickness attains 250 nr in the spatial domain. 
The sequence encompasses the Bajocian (Dl) and Bathonian limestones (D3 
and D4) representing the main aquifer formations of the Dogger, a clastic 
dominated interval (’’Marries de Longwy”, D2) separating the two. The Dog¬ 
ger sequence is topped with a thin oolithic limestone from Lower Callovian 
(’’Dalle Nacree”, Cl), implemented as a 15 m thick layer in the model. The 
latter marks the transition with the thick, highly impermeable, claystone for¬ 
mation of Callovo-Oxfordian age (C2) of which the thickness reaches 150 m 
in the model. In the numerical simulations, a target zone (TZ) located in the 
middle of layer C2 (Figure [l]) represents the location for the computation of 
the output quantity of interest. 

The low-permeability COX layer is overlaid by a limestone sequence of 
the Oxfordian age. The latter is incorporated as a 260 m thick formation sub¬ 
divided into 6 hydrogeological entities. A relatively confining layer from the 
Upper Argovian (C3ab) rests directly on the COX and is followed by perme¬ 
able formations of the Rauracian-Sequanian sequence (Lla to L2c). A thick 
interval of marls and argillaceous limestones from Kimmeridgian age (Kl- 
K2) covers the whole and is implemented as a 160 m thick low-permeability 
layer. The top layer is a 120 m thick confining formation attributed to the 
Tithonian (K3). The latter outcrops in the vicinity of Bure. 

2.2 Governing equations and model outputs 

In the numerical simulations the flow is governed by the steady-state equation 

V ■ q = 0, (1) 

where q = — K ViL, is the Darcian flux vector [L T -1 ]. K is the tensor of hy¬ 
draulic conductivity [L T -1 ] and H is the hydraulic head [L], The anisotropy 
Ak in the components of the tensor of hydraulic conductivity is defined as 
the ratio between the hydraulic conductivities in the two principal directions: 
A k = K z /K x . 

Here, it is assumed that K has orthotropic properties. Considering a 
hydraulic conductivity tensor K p of which the components are mapped into 
the Cartesian system and given along their principal direction, X p , the tensor 
K in the global Cartesian space is retrieved by means of the rotation matrix 
R with the expression 

K = R t K p R. (2) 

For the two-dimensional problem considered, the rotation matrix R is defined 
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Figure 1: Geometry and geological layers with the localization of the target 
zone (vertical exaggeration: 20). 


in terms of the Euler angle 9 (in degree) as 


R 


cos 9 sin 9 
- sin 9 cos 9 


( 3 ) 


In the present study, steady-state flow simulations are carried out to¬ 
gether with the computation of the lifetime expectancy probability density 
function (PDF) at any point x in the domain. Under stationary conditions 
(i.e. steady-state flow), the lifetime expectancy PDF addresses the proba¬ 
bility distribution of the time required for a solute, taken at any position 
x, to leave the domain. In its formulation, the lifetime expectancy PDF 
assimilates the forward advective-diffusive transport equation (ADE) to the 
Fokker-Planck (forward Kolmogorov) equation measuring the random motion 


of solute particles (Uffink, 1989). For more details on the computation of the 


lifetime expectancy PDF, the reader is referred to Cornaton and Per rochet 
( 2006a|b ) and Kazemi et al. (2006). 

Based on the ADE, the lifetime expectancy PDF is computed using the 
backward transport equation requiring reversed flow directions (q := —q) as 
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well as adjusted downstream boundary conditions. The lifetime expectancy 
PDF g£;(x,t) at any point x in the domain is then governed by 

= v ■ (q g E + D V 9b ), (4) 

where 0 is the effective porosity [-] and where D is the dispersion tensor 

0D = (a L - + aT H^ll 1 + $ D m I, (5) 

where I is the identity matrix, D m is the coefficient of molecular diffu¬ 
sion [L 2 T -1 ], aL and ax are the longitudinal and transverse components 
of the macro-dispersion tensor [L] respectively. In the present study, the 
anisotropy in the macro-dispersion tensor is determined with the coefficient 
A a = ax/ o>l- 

The straightforward computation of the first moment of the lifetime ex¬ 
pectancy PDF is the so-called mean lifetime expectancy (MLE) E(x,z) [T] 
at any position x, governed by 


V • (q E + D VE) = 0, 


( 6 ) 


where it can be seen that the porosity 0 [-] acts as the sink term in the aging 
process. 

The TZ comprises a set of 1,947 nodes in layer C2, covered by a rectangle 
of which the lateral and vertical extensions are x = [18,440; 21,680], z = 
[425; 435] (Figure [l]). In the present study, the arithmetic mean of E(x,z ) 
calculated at each of these 1,947 nodes stands for the output response of 
interest and is used in the subsequent analysis. It can be seen as the average 
time for a solute originating from the TZ to reach any outlet of the domain 
of the model. 

The finite element simulator GroundWater (Cornaton 2007) was em¬ 


ployed to solve Eq. Q-Q using the finite element technique. A single run 
of steady-state flow and MLE computation takes about 120 seconds using a 
parallel solver with 6 CPU. 

The reader should note that the use of a 2D vertical model to solve for 
the hydro-dispersive processes cannot capture correctly the real behavior of 
the Paris Basin subsurface because, apart from being a simplified model, it 
omits the lateral flow and dispersion along the third dimension. This has the 
effect of underestimating the magnitude of the modeled processes[Kerrou and 
Renard (2010). It is however assumed that this bias is equivalent for all the 


layers considered and thus, the interpretation of the SA results obtained with 
the 2D cross-section may be generalized to a synthetic 3D case employing 
the same settings. 













2.3 Flow boundary conditions 


The fully saturated model considers stationary flow conditions in a confined 
aquifer which are implemented as Dirichlct type flow boundary conditions 
(BCs). These flow BCs are imposed on nodes located at the top of the 
model domain as well as at both lateral limits of the two limestone sequences 
(Figure [2]). 


Regional piezometric maps based on field measures (Linard et al. 2011) 


were used to constrain the hydraulic gradients in both carbonate sequences. 
The flow BCs imposed on the lateral boundaries of the two limestone se¬ 
quences derive from a 25 km transect starting from the Gondrecourt trough 
and extending in a North-West direction, the main regional flow direction. 
The hydraulic gradient set on top of the domain corresponds to the average 
topographic gradient of the region covered by the transect. 

Under these conditions, the general groundwater flow direction is oriented 
from right to left. The proportions of the total outflowing rates are approx¬ 
imately 2%, 60% and 38% for the top of the domain, the Oxfordian and 
the Dogger discharge boundaries, respectively. In layer C2, the groundwater 
flows downward in the very right part of the domain and then upward in the 
remainder, with a hydraulic gradient inversion in the vicinity of the TZ (see 
Figure [ 2 ]). As a summary, the flow BCs are gathered in Table [I] 


Table 1: Flow boundary conditions. 


Boundary 

Position 

Hydraulic head 

right Oxfordian 

x = 25000, 2 = [500, 760] 

H = 305 m 

left Oxfordian 

x = 0, z = [500, 760] 

H = 230 m 

right Dogger 

x = 25000, z = [110, 360] 

H = 295 m 

left Dogger 

x = 0, z = [110, 360] 

H = 275 m 

top of domain 

x = [0, 25000], z = 1040 

H = 225 + 85:r/25000 

elsewhere 


no flow 


To account for uncertainties in the flow BCs, the hydraulic gradients in 
the two limestone sequences and on the top of the domain are considered as 
uncertain input factors included in the following SA (Section [4]). A change 
in the hydraulic gradients may shift the position of the vertical groundwater 
flux inversion in layer C2, and thus the MLE calculated at the TZ. 

2.4 Hydraulic conductivity and porosity values 

Many studies have undertaken the inventory of hydraulic conductivity (A') 
and porosity (</>) values in the various geological formations of the Paris Basin. 
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H = 225 m H = (225 + 85 l / 25000) m H = 310m 
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Figure 2: Flow boundary conditions and head contours (vertical exaggera¬ 
tion: 20). 


For a large number of wells and boreholes within a wide area around the 
experimental site of Bure, laboratory and field measurements were conducted 



1 2011; 

Fourre et al. 

2011; Delay and Distinguin 

2004 


However, very few {K, 0} datasets are available for the four low-permeability 
formations implemented in the present model (he. K3, K1-K2, C2 and T). 


Hence, data extracted from the literature (Cosenza et ah, 2002 Delay et ah 


2007b!, 

2006; 

Enssle et ah, 

2011 


Mazurek et al. 


and employed in previous studies (Contoux et al. 


2011 Vinsot et al. 


2013 de Hoyos et al. 


2011 ), 


2012 


'or the 


|Goncalves et al. , 2004a|b ), were used to define the uncertainty ranges 
{A', 0} sets in these layers. 

In the geological formations of the Oxfordian and Dogger sequences, large 
variabilities of the {A', 0} couples are noticed with the presence of dependen¬ 
cies ( e.g. low K and low 0 values are correlated). However, in order to 
simplify the conceptual approach in a first stage, a perfect dependence be¬ 
tween log lo (A0 and 0 is defined by making use of a mathematical function 
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approximating the relationship between these two in each layer. In the se¬ 
quel, both parameters are referred to as a whole under the name petrofacies. 
This approach reduces the computational burden of the subsequent SA by 
avoiding the use of correlation functions between the two uncertain factors. 
The hydraulic conductivity values deriving from the {K , 0} distributions in 
each layer are attributed to the longitudinal component of the hydraulic con¬ 
ductivity tensor, K x . For each layer, the estimated value of K x is retrieved 
through a relationship: log 10 (A' x ) = /(0) (Deman et al. 2015). 

Although no explicit information is available on the following, the geo¬ 
logical formations are believed to feature anisotropic hydraulic conductivity 
tensors K. i.e. anisotropy in the two principal components of the tensor (K x 
and K z ) defined as the ratio Ax = K z /K x . In the nominal case, Ax = 0.1 
is assumed for every layer in the model. 

Preferential flow directions are supposedly taking place within each indi¬ 
vidual layer. For each layer, the Euler angle 6 defines the orientation of the 
hydraulic conductivity tensor ~K P in the Cartesian space (see Eq. (§-(§)• In 
the nominal case, 9 = 0 degree is assumed in every layer, which corresponds 
to the two principal components of the hydraulic conductivity tensors 
being oriented along the x and z axes. The orientation of the groundwater 
flux q in the model is principally due to the static hydraulic gradients VH 
resulting from flow BCs implemented on the edges of the domain. Note how¬ 
ever that the Euler angle 6 may locally change the orientation of q in a given 
layer and thus, drive the groundwater into adjacent layers where magnitudes 
might be different. This phenomenon may have a significant effect on the 
MLE calculated from the TZ, which is explored in Section |4j 

Table [2] summarizes the nominal values for 0 and the corresponding K x 
in each of the 15 hydrogeological layers comprised in the model. The values 
for 0 correspond to the mean value (or the median value of the CDF) of the 
distribution in each layer, whereas the values for K x derive from approxi¬ 
mation functions. As a reminder, the present study assumes homogeneous 
parameters in every layer. Although this feature is unrealistic, it is recalled 
that the purpose of this study is to bring insights into the global effect of 
equivalent advective-dispersive parameters of the multi-layered hydrogeolog¬ 
ical model and to provide recommendations for a similar application on a 
high-resolution integrated hydrogeological model of the Paris Basin. 

2.5 Dispersion parameters 

The mean lifetime expectancy formulation (Eq. (|6])) is an advective-dispersive 
solute transport equation, where the longitudinal and transverse components 
of the macro-dispersion tensor (oll and ar respectively) control the particles 


(Deman et al. 2015 
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Table 2: Nominal values for the porosity (</>) and the longitudinal hydraulic 
conductivity (K x ) in the 15 hydrogeological layers. 


Layer K x [m/s] (f> 


K3 

9.01E—09 

0.0100 

K1-K2 

4.53E-09 

0.1150 

L2c 

1.10E—06 

0.1389 

L2b 

3.46E-07 

0.1110 

L2a 

1.62E—07 

0.1139 

Lib 

1.49E-05 

0.1604 

Lla 

1.17E-06 

0.1549 

C3ab 

4.59E-08 

0.0984 

C2 

1.99E-13 

0.1580 

Cl 

1.89E-06 

0.0470 

D4 

1.65E-05 

0.0905 

D3 

1.76E-06 

0.1016 

D2 

2.62E-07 

0.0623 

D1 

3.23E-06 

0.0688 

T 

1.95E-12 

0.0810 


dispersion. These two uncertain factors depend particularly on the rock 
type, on the tortuosity of the porous media and also, on the scale considered. 
Homogeneous values of c\j J and olt are set within each layer, with the values 
otL = 15 m and A a = olt/oll = 0.1 considered in the entire domain of the 
model in the nominal case. 

As mentioned previously, no decay or adsorption effects are accounted 
in the computation of the MLE. The coefficient of molecular diffusion cor¬ 
responds to the theoretical self-diffusion coefficient for the water molecule, 
D m = 2.3 10“ 9 m 2 /s. 

3 Polynomial chaos expansions for sensitivity 
analysis 

Let us denote by M. the computational model describing the behavior of 
the considered physical system. Let X = {Xi, ... , Xm } denote the M- 
dimensional random input vector with joint PDF fx (x) and marginal PDFs 
fxi( x i)i i — 1, • • • , M. Due to the input uncertainties represented by X, 
the model response of interest becomes random. The computational model 
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is thus considered as the map 

X G V x C M m i —> Y = M(X) e M, (7) 

where is the support of X. In the description of the theoretical frame¬ 
work hereafter, we assume that the components of X are independent , which 
is the case for the model in the present study. 

As explained in the Introduction, the aim of global SA is to identity ran¬ 
dom input variables and combinations thereof with significant contributions 
to the variability of Y as described by its variance. A concise description of 
the employed method of PCE-based Sobol’ sensitivity indices is given in the 



3.1 Sobol’ indices 


Assuming that the function JYl(X) is square-integrable with respect to the 
probability measure associated with fx(x ), the Sobol’ decomposition of Y = 
A4 (X) in summands of increasing dimension is given by Sobol’ (1993) 

M 

M(X) = M 0 + J2Mi{X 1 )+ M ij {X i ,X j ) + ... + M 12 ... M (X) (8) 

i =1 1 <i<j<M 

or equivalently, by 


M(X)=M 0 + Y, M n( x u), 


(9) 


where Afo is the mean value of Y, u = {A, ... , C {1, ... , M} are index 
sets and X u denotes a subvector of X containing only those components 
of which the indices belong to u. The number of summands in the above 
equations is 2 M — 1. 

The Sobol’ decomposition is unique under the condition 



M u (x u )fx k (x k )dx k = 0, 


if k E u, 


( 10 ) 


where T>x k and fx k (x k ) respectively denote the support and marginal PDF 
of X k . Eq. (10) leads to the orthogonality property 

E [M. U (X U )M. V (X V )\ = 0, if u^v. (11) 
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The uniqueness and orthogonality properties allow decomposition of the vari¬ 
ance D of Y as 

D = Vm[M(X)] = Y j D u , (12) 

where D u denotes the partial variance 

D u = Var [M U (X U )\ = E [M 2 U (X U )] . (13) 

The Sobol’ index S u is defined as 

S u = D u /D, (14) 

and describes the amount of the total variance that is due to the interac¬ 
tion between the uncertain input parameters comprising X u . By definition, 
= 1. First-order indices, Si, describe the influence of each pa¬ 
rameter X, considered separately, also called main effects. Second-order in¬ 
dices, Sij, describe the influence from the interaction between the parameters 
{X i 7 Xj}. Higher-order indices describe influences from interactions between 
larger sets of parameters. 

The total sensitivity indices, Sf, represent the total effect of an input 
parameter X t , accounting for its main effect and all interactions with other 
parameters. They are derived from the sum of all partial sensitivity indices 
S u that involve parameter X t , i.e. 

Sj = J2 D ™/ D i = (15) 

ii 

It follows that ST = 1 — S~i, where S~i is the sum of all S u with u not 
including i. 

Evaluation of the Sobol’ indices by Monte Carlo simulation is based on a 
recursive relationship that requires computation of 2 M Monte Carlo integrals 
involving XI(X) in order to obtain the full set. Clearly, this is not affordable 
when the computational model is a time-consuming algorithmic sequence. In 
typical applications, first-order and total indices are computed. We note that 
the first-order indices are equivalent to the first-order indices obtained by the 
FAST method, which may provide a more efficient computation. It will be 
shown next that when a PCE of the quantity of interest is available, the full 
set of Sobol’ indices can be obtained analytically at almost no additional 
computational cost. 
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3.2 Polynomial chaos expansions 


3.2.1 Computation of polynomial chaos expansions 

A PCE approximation of Y = JA{X) in Eq. ([T]) has the form Xiu and 


Karniadakis (2002) 


Y = M PCE (X) = ^ VMX), (16) 

aeA 


where A is a set of multi-indices ct = (cq, ... ,chm), G A} is a set 

of multivariate polynomials that are orthonormal with respect to fx(x) and 
{y a ,ct G A\ is the corresponding set of polynomial coefficients. 

The multivariate polynomials that comprise the PCE basis are obtained 
by tensorization of appropriate univariate polynomials, i.e. 

M 

S'„(X)=n^(V i ), (17) 

i=l 


where t/Y (Xj) is a polynomial of degree cq in the i-th input variable be¬ 
longing to a family of polynomials that are orthonormal with respect to 
fxAxi). For standard distributions, the associated family of orthonormal 
polynomials is well-known e.g. a standard normal variable is associated with 
the family of Hermite polynomials, whereas a uniform variable over [—1,1] 
is associated with the family of Legendre polynomials. A general case can 
be treated through an isoprobabilistic transformation of the input random 
vector X to a basic random vector e.g. a vector with independent standard 
normal components or independent uniform components over [—1,1]. The 


set of multi-indices A in Eq. (16) is determined by an appropriate truncation 
scheme. In the present study, a hyperbolic truncation scheme is employed, 
which corresponds to selecting all multi-indices that satisfy 



S P- 


(18) 


for appropriate 0 < q < 1 and p G Id Blatman and Sudret (2010b). When 


q = 1, polynomials of maximum total degree p are retained, whereas use of 
a lower q limits the number of basis terms that include interactions between 
two or more variables. Optimal values of these parameters may be selected 
in terms of error estimates e.g. by using cross validation techniques. 

Once the basis has been specified, the set of coefficients y = {y ai ct G A} 
may be computed by minimizing the mean-square error of the approximation 
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over a set of N realizations of the input vector, £ = (ad 1 ), ... , ad 7 ^}, called 
experimental design. Efficient solution schemes are obtained by considering 
the regularized problem 

M(x^) - ]T v a # a (x®)) + CIMI?, (19) 

aeA ' 


N 

y = arg min 

ueK card - A 


where ||t;|| 


E card A i 
.7—1 Inl¬ 


and C is a non-negative constant. A nice feature 


of the above regularized problem is that it provides a sparse meta-model by 
disregarding insignificant terms from the set of predictors. Higher values of 
C lead to sparser meta-models, while its optimal value is typically identified 
as the one leading to the minimum error estimated with e.g. cross validation 
techniques Tibshirani (1996). In the present application, we solve Eq. (19) 


using the hybrid Least Angle Regression (LAR) method as originally proposed 


in 

Blatman and Sudret 

(2011 

). Hybrid LAR employs the LAR algorithm 

Efron et al. 

(2004 

) to select the best set of predictors and subsequently, 


estimates the coefficients with standard least-squares minimization. 


3.2.2 Error estimates 

A good measure of the accuracy of PCE is the mean-square error of the 


residual, Errc = E 


Y - Y 


, called generalization error. In practice, 


this could be estimated by Monte Carlo simulation using a sufficiently large 
set of n va i realizations of the input vector, A va i = {aq, ... ,a? nval }, called 
validation set. The estimate of the generalization error is given by 


^val / 

Err g =- Y' ( M{xi) - V' y a ^ a (xi 

n v al “ V 

1=1 CX.£A 


( 20 ) 


The relative generalization error, erre, is estimated by normalizing Errc 
with the empirical variance of 3^ V ai = {Al(aq), ... ,J\4(x Uva J}, denoting the 
set of model evaluations at the validation set. 

However, PCE are typically used as surrogate models in cases when evalu¬ 
ating a large number of model responses is not affordable. It is thus desirable 
to get an estimate of Errc using only the information obtained from the ex¬ 
perimental design. One such measure is the Leave-One-Out (LOO) error 
Allen (1971). The idea of the LOO cross-validation is to set apart one point 
of the experimental design, say an*q and use the remaining points to build 
the PCE, denoted A1 PCE ^. The LOO error is obtained after alternating over 
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all points of the experimental design, i.e. 


1 N 

Err LOO = ^ £( M(x (<) ) - M pce \*(xV)) 2 . ( 21 ) 

i=1 

Although the above dehnition outlines a computationally demanding pro¬ 
cedure, algebraic manipulations allow evaluation of the LOO error from a 
single PCE based on the full experimental design. Let us denote by hj the i- 
tli diagonal term of matrix where th = {^ij = \Etj(a;W), i = 

1, ... , N; j — 1, ... , card *4.}. Then, the LOO error can be computed as 


Err loo = 


1 A /— ,M PCE (ccW) 

V 1 -hi 

1=1 x 


2 


( 22 ) 


The relative LOO error, errioo, is obtained by normalizing ErrLoo with 
the empirical variance of y = {JA(x^), ... , Ad(a^A*)}, denoting the set of 
model evaluations at the experimental design. Because this error estimate 
may be too optimistic, a corrected estimate is herein employed, given by 


Chapellc et al. (2002) 


err loo = err LO o 


cardAA 

) 


-l 

(1 + tr(('I' T 'I')- 1 )) . 


(23) 


This corrected LOO error is a good compromise between fair error estimation 
and affordable computational cost. 


3.3 Sobol’ indices from polynomical chaos expansions 

Let us consider the PCE Y = Ad PCE (X) of the quantity of interest Y = 
A4(X). It is straightforward to obtain the Sobol’ decomposition of Y in an 
analytical form by observing that the summands M„ CE (X„) in Eq. (J9|) can 
be written as 

Ml CE (X u ) = Y, 3/«*«, (24) 

ql£>A. u 

where A u denotes the set of multi-indices that depend only on u , i.e. 

A u = (a £ 4 : Oji ^ 0 if and only if k £ it}. (25) 

Clearly, (J A u = A. Consequently, due to the orthonormality of the PCE 
basis, the partial variance D u is estimated as 

D u = Var [Ml CE (X u )} = £ yj, (26) 


17 







and the total variance is estimated as 

D = Var [M PCE (X)j = ^ y a 2 . (27) 

«eyl\{o} 

Accordingly, the Sobol’ indices of any order can be obtained by a mere com¬ 
bination of the squares of the PCE coefficients. For instance, the PCE-based 
first- and second-order Sobol’ indices are respectively given by 

Sj= Vo?/D, Ai = {a E A: ai> 0, aq/; = 0} (28) 

aeAi 

and 

Sij ^2 y<x 2 /D, Aij = {a G A : a u aq > 0, a k ^ id = 0}, (29) 

ctCAij 

whereas the total Sobol’ indices are given by 

Sf = ^2 A? = {a <E A: ai> 0}. (30) 

aeAf 

ft is evident that once a PCE representation of Y = Ai(X) is available, 
the complete list of Sobol’ indices can be obtained at a nearly costless post¬ 
processing of the PCE coefficients requiring only elementary mathematical 
operations. 

4 Results and discussion 

Figure [3] provides an overview of the distribution of the MLE throughout the 
entire model domain in the nominal case. In this case, the parameters are 
distributed homogeneously in each of the 15 layers; K x and (j) take on the 
values given in Table [2] for each layer, whereas for all layers the anisotropy 
ratio is Ak = 0.1, the Euler angle is 6 = 0 degree, the longitudinal component 
of the macro-dispersion tensor is = 15 nr and the anisotropy ratio is 
A a = 0.1. The hydraulic gradients follow the boundary conditions settings 
described in Section 12.31 

Because of its highly confining properties, the middle layer (C2) presents 
values of MLE > 40,000 years. On average, it takes approximately 75,000 
years for a solute departing from the TZ to reach any outlet of the domain. 
Much lower MLE values are found in the two aquifer sequences, with the 
Oxfordian displaying slightly smaller values. The effect of conductive layers 
is clearly distinguishable as fringes of low MLE values stretch in layers D4, 
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Figure 3: Spatial distribution of mean lifetime expectancy in the reference 
case (vertical exaggeration: 10). 


Lla and Lib in particular. As a result of the low permeability in the top two 
layers (K3 and K1-K2) and the bottom layer (T), water molecules can take 
more than a 100,000 years to flow through the domain. 

In the following, we compute the PCE-based Sobol’ indices for the MLE 
at the TZ by implementing the theory presented in Section [3j We con¬ 
sider a high-dimensional random input encompassing the entirety of hydro- 
dispersive parameters described in Section [2] as well as the flow boundary 
conditions. Note that the term porosity , 0, is construed in the discussion 
of the Sobol’ indices. Since the values of the hydraulic conductivities are 
retrieved through approximation functions, the estimation of the sensitivity 
for the 0 variables is implicitly associated to that of the K x variables in the 
respective layers. This is singularly important when interpreting the Sobol’ 
indices for aquifer formations, where the hydraulic conductivity governs the 
ageing process (Section 4.4). Computations of the PCE and Sobol’ indices 


are performed with the uncertainty quantification software UQLab Marelli 


and Sudret (2014); Marelli et al. (2015). 


4.1 Modeling of input uncertainty 

The uncertain input in each of the 15 layers of the hydrogeological model 
comprises the following parameters governing the advective-dispersive pro¬ 
cesses: the petrofacies, P; the anisotropy in the components of the hydraulic 
conductivity tensor, Ax] the Euler angle of the hydraulic conductivity ten¬ 
sor, 6 ; the longitudinal component of the dispersivity tensor, a L ; and the 
anisotropy in the longitudinal and vertical components of the dispersivity 
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tensor, A a . Furthermore, the uncertain input includes the hydraulic gradi¬ 
ents, ViL, at three zones: the Dogger sequence, the Oxfordian sequence and 
the top of the model domain, respectively denoted as zone 1, zone 2 and 
zone 3. 

As stated before, a deterministic relationship log 10 (/i x ) = /(0) is assumed 
for each layer, i.e. the uncertainty regarding the petrofacies, P, of a layer 
is treated through the porosity, 0. The available sparse data, characterizing 
specific locations in the subsurface of the Paris Basin, provide insufficient 
information for a comprehensive definition of the CDF of the 0 parameters. 
Thus, in order to conduct a conservative study and avoid introducing any 
bias in the analysis, the uncertain porosities are modeled as uniform random 
variables. The porosity ranges are bounded by the values 0( mm ) and 0( m< “) 
listed in Table [3] together with the respective coefficients of variation (CoV); 
these are shown graphically along the model cross-section in Figure [4} The 
bounds 0( mm ) and 0( maa: ) represent the 1st and 9th deciles of the correspond¬ 
ing CDF derived from porosity values measured in each geological layer. This 
approach is justified by the presence of local extreme measures that cannot 
be representative for the whole layer. Bounds for the K x parameters are also 
provided in Table [3j consistently with the log lo (A0) = /(0) approximation 
functions, and presented graphically in Figure [5] 

Table 3: Ranges of porosity, 0, and respective CoV and hydraulic conduc¬ 
tivity values, K x , in the 15 geological layers. 


Layer 0( mm ) 


max ) 

-] CoV K [ ™ n) 

m/s 

^-(max) 

m/s 


K3 

0.0840 

0.1160 

0.0924 

3.3734e-10 

2.4078e-07 

K1-K2 

0.0870 

0.1430 

0.1406 

9.8116e-ll 

2.0928e-07 

L2c 

0.1019 

0.1759 

0.1538 

3.6186e-08 

2.6212e-06 

L2b 

0.0645 

0.1574 

0.2417 

8.7318e-10 

6.3950e-06 

L2a 

0.0651 

0.1627 

0.2474 

4.7005e-10 

9.9336e-06 

Lib 

0.1375 

0.1833 

0.0824 

3.4324e-09 

2.8913e-04 

Lla 

0.0991 

0.2107 

0.2080 

3.1165e-08 

2.1523e-06 

C3ab 

0.0747 

0.1221 

0.1391 

7.8488e-09 

1.2945e-06 

C2 

0.1284 

0.1876 

0.1082 

5.0349e-14 

6.2570e-13 

Cl 

0.0142 

0.0799 

0.4031 

1.8184e-07 

1.6195e-05 

D4 

0.0237 

0.1573 

0.4262 

1.6408e-07 

3.1521e-03 

D3 

0.0237 

0.1795 

0.4427 

1.7470e-07 

4.3539e-06 

D2 

0.0185 

0.1061 

0.4059 

6.6071e-08 

1.7049e-06 

D1 

0.0191 

0.1186 

0.4172 

6.2552e-08 

1.8425e-05 

T 

0.0696 

0.0925 

0.0816 

1.2325e-13 

8.1328e-12 
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Figure 4: Porosity ranges along the model cross-section. 


Due to the lack of explicit information, each of the parameters Ak , 0, oil 
and A a is identically distributed in the different layers. In particular, the 
anisotropy ratios Ak and A a both follow a uniform distribution in [0.01, 1], 
the Euler angle 6 is uniformly distributed in [-30, 30] (in degrees) and the 
parameter oll is uniformly distributed in [5, 25] (in meters). The static hy¬ 
draulic gradients are also uniformly distributed in the ranges given in Table 
[4} These were obtained by applying a perturbation of 20% of the nominal 
hydraulic gradient within each of the three zones. Although hypothetical, 
these conservative uncertainty ranges were purposely chosen to provide in¬ 
sights into the behavior of the multi-layered model. 

Table 4: Ranges of hydraulic gradient at the three zones of interest. 


Zone number 

V [{(.min) 

y ][{max) 

1 

0.00064 

0.00096 

2 

0.00240 

0.00360 

3 

0.00272 

0.00408 


In total, the model subsumes M = 78 independent input random variables 
upon the MLE of water molecules outflowing from the TZ in the middle of 
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Figure 5: Hydraulic conductivity ranges along the model cross-section. 


layer C2. 


4.2 Polynomial chaos expansions of the model response 

In the sequel, we build PCE of the model response in terms of the 78 input 
random variables. To this end, we employ Latin hypercube sampling (LHS) 
to draw the experimental design McKay et al. (1979). LHS is a popular 


technique for obtaining random experimental designs ensuring uniformity of 
each sample in the domains of the margin input variables {Ad, ... ,Xm}- 
The data available herein for building the PCE consist of the MLE values for 
(i) an experimental design of size N = 2, 000, drawn with LHS and denoted 
£, and (ii) an enrichment of £ of equal size N' = 2, 000, denoted £'. The 
enrichment is built according to the nested-LHS approach so that the joint 


set {£,£'} is approximately a LHS experimental design as well Wang (2003); 


Blatman and Sudret (2010c). Histograms of the model response for the two 


sets of input vectors are shown in Figure [6} Positively skewed distributions 
are observed for both output sets with the modes situated at MLE ~ 85,000 
years. 

We develop PCE based on £ and on the joint set {£, £'} and assess their 
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Figure 6: Histograms of mean lifetime expectancy values calculated with sets 
S (left) and S' (right). 


comparative accuracy. For the set S, we consider the MLE response in both 
the original and the logarithmic scales; in this case, the enrichment S' serves 
as a validation set for computing the generalization error (see Section 3.2.2). 
For all PCE, the candidate basis is determined using a hyperbolic truncation 
scheme with q = 0.5 (see Eq. (18)). Sparse PCE are developed by varying 
the maximum degree p from 1 to 15; the meta-model of optimal degree is 
selected as the one yielding the smallest corrected LOO error (see Eq. (23)). 

The first PCE, denoted A, is built using S as the experimental design 
and considering the MLE in the original scale. The optimal degree is p = 8 
and the corresponding corrected LOO error is err* LOO = 0.0565. The sparse 
PCE includes 185 basis elements, whereas the total number of basis elements 
for p = 8 and q = 0.5 is 18,643; for q — 1, the size of the candidate basis 
would be 5.3 x 10 10 . The index of sparsity, defined as the number of basis 
elements in the sparse expansion divided by the size of a full basis for the 
same p and q, is 185/18,643 = 9.9 x 10 -3 . The small value of this index 
herein indicates the interest in developing sparse PCE for analyses of high¬ 
dimensional models. The sparse basis consists of polynomials in 68 out of the 
78 total input parameters, meaning that the output does not depend at all on 
the values of the 10 excluded parameters. Note that 3 out of the 10 excluded 
parameters are properties of layer T. The estimate of the generalization error 
(evaluated with S') is erre = 0.0759. The left graph of Figure [7] compares the 
values of the meta-model, Y, with the respective values of the exact model, 
Y, at the input samples of the experimental design, S. A similar comparison 
but for the validation set, S ', is shown in the right graph of the same figure. 

A second PCE, denoted B, is built by using again S as the experimental 
design, but employing a logarithmic transformation of the MLE. The opti- 
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errc=0.07591 




Figure 7: Comparison of PCE A with the actual model response at the 
experimental design, £, (left) and at the validation set, S' (right). 


mal sparse PCE is obtained for p = 8 (same degree as for PCE A) with 
corresponding corrected LOO error err* LOO = 0.0287. It comprises 163 ba¬ 
sis elements and thus, has an index of sparsity 163/18,643 = 8.7 x 10~ 3 . 
The sparse basis consists of polynomials in 65 out of the 78 total input pa¬ 
rameters, with 6 out of the 13 excluded parameters representing dispersivity 
anisotropy ratios (A a ). The resulting generalization error for the MLE re¬ 
sponse in the original scale (considering the exponential transformation of 
the obtained PCE) is estimated as errc = 0.0452. Note that both err* LOO 
and errc are lower than the respective error estimates for PCE A. The left 
and right graphs of Figure [8] compare the exponential transformation of PCE 
B with the exact model response at £ and £', respectively. 

Finally, we use as experimental design the joint set {£,£'} consisting of 
N + N' — A, 000 points. The optimal PCE is obtained for p — 10 and the 
corresponding corrected LOO error is err* LOO = 0.0384. The sparse PCE 
comprises 312 basis elements and has an index of sparsity 312/106,887 ~ 
2.9 x 10” 3 . The only two parameters excluded from the sparse basis are Aj^ b 
and «£. The comparison between this meta-model, denoted C, and the exact 
model response at the input samples of the experimental design, {5,5'}, is 
shown in Figure [9j 

Assessing the relative accuracies of the three meta-models, we note that 
all have (corrected) LOO errors of the same order of magnitude, with the 
smallest error corresponding to PCE B. Because it is of interest to limit the 
number of costly evaluations of the exact hydrogeological model, an exper¬ 
imental design comprising 2, 000 points is deemed most appropriate. We 
therefore conduct SA for the MLE response by post-processing the coeffi¬ 
cients of PCE A and B and compare the results. We underline that the 
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Figure 8: Comparison of exponential transformation of PCE B with the ac¬ 
tual model response at the experimental design, £, (left) and at the validation 
set, S' (right). 



Figure 9: Comparison of PCE C with the actual model response at the 
experimental design, {£,£'}. 


Sobol’ indices are not invariant under a logarithmic transformation of the re¬ 
sponse quantity of interest (see e.g. Borgonovo et al. (2014)). Nevertheless, 


we expect that the Sobol’ indices for the non-transformed and the logarith¬ 
mic MLE will exhibit similar trends, i.e. a parameter that is important for 
the variance of the MLE will also be important for the variance of its loga¬ 
rithm. Moreover, it is of interest to perform SA with PCE B, because the 
logarithmic transformation represents herein a meaningful quantity from a 
physical viewpoint considering the wide variation of MLE. 

To further elaborate on the accuracy of the two meta-models employed 
in the subsequent SA, we compare the ED-based mean, /i, and standard 
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deviation, a, of the response with the respective quantities, /2 and a, obtained 
from the PCE coefficients. For the response in the original scale (PCE A), 
we have {/i = 113,216; a = 53,220} and {/i = 113,410; a = 50,735}, 
corresponding to errors (e^ = 0.17 %; = —4.67 %}. For the logarithmic 
response (PCE B), we have {/i = 11.545; cr = 0.4180} and {jj, = 11.543; a = 
0.4073}, leading to errors (e^ = —0.01 %; e$ = —2.56 %}. Overall, PCE B 
is more accurate. 


4.3 Sensitivity analysis 


Figure [10] shows bar-plots of the total and first-order Sobol’ indices for PCE 
A. The ten largest indices are presented in descending order. The superscripts 
on the parameter symbols on the horizontal axes denote layer names or zone 
numbers. The figures indicate the same ranking of the Eve most important 
parameters in terms of both the total and first-order indices. These are the 
porosities of layers D4, C3ab, Lib, Lla and Cl in order of importance, with 
the porosity of layer D4 being dominant. Employing the criterion Sj < 0.01 
to sort out unimportant variables, the porosities of the aforementioned five 
layers comprise the only important parameters. Therefore, the screening 
allows one to consider 73 out of 78 parameters as unimportant, meaning 
that they could be given a deterministic value without affecting essentially 
the predicted MLE. 

We note that all five layers with porosities identified as important are 
located close to the host layer C2. D4 is the thickest among those and has 
the highest hydraulic conductivity. Although Cl is adjacent to the host 
layer C2, it is associated with smaller total and first-order indices than other 
neighboring layers, which may be attributed to its small thickness. 

The largest second-order effects (not shown) comprise interactions be¬ 
tween the five porosities classified above as important and involve one of the 
layers D4 or Lib. The second-order effects explain 10.9% of the total vari¬ 
ance, whereas contributions from higher-order effects are practically zero. 

In the sequel, we demonstrate that the sensitivity pattern in this high¬ 
dimensional problem can be detected by only a few runs of the model. To 
this end, we randomly extract 200 points out of the experimental design E 
and follow the same procedure as above to compute the PCE-based sensi¬ 
tivity indices for the MLE in the original scale. Conducting 100 repetitions, 
we obtain statistics for the ten largest total Sobol’ indices identified in the 
analysis with PCE A (see upper graph of Figure 10). These are depicted 


in the boxplot in Figure 11 in which the central mark of a box indicates 


the median and its edges indicate the 25th and 75th percentiles. We observe 
that the median Sobol’ indices from SA performed with experimental designs 
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Figure 10: Sobol’ indices using PCE A. 


comprising only 200 points are fairly close to the Sobol’ indices obtained with 
the experimental design £ comprising 2,000 points. Moreover, the dispersion 
is relatively small, i.e. ±0.03 at most, which is more than sufficient for a 
preliminary screening of the important variables. These results emphasize 
the robustness and remarkable efficiency of our approach in dealing with 
high-dimensional SA problems. 

Let us now compare the above results with respective results obtained 
by considering the logarithmic MLE. In Figure [I2j we show bar-plots of the 
ten highest total and first-order Sobol’ indices in descending order for PCE 
B. Obviously, the indices obtained in the two cases, i.e. considering the 
non-transformed and the logarithmic MLE, follow similar trends. SA with 
PCE B identifies the same five variables as important (Sf > 0.01) in the 
same order as the analysis with PCE A, with the porosity of layer D4 being 
dominant. Furthermore, the Eve important variables are also the ones with 
the highest first-order indices following the same ranking. Again, the largest 
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Figure 11: Boxplot of total Sobol’ indices for the MLE in the original scale 
using experimental designs with 200 points. 


second-order effects involve the porosities of one of the layers D4 or Lib, while 
their sum explains 4.2% of the total variance. Thus, second-order effects are 
slightly less significant for the logarithmic response than for the response 
in the original scale. Contributions of higher than second-order effects are 
practically zero. 

To gain further insight into the effects of the important variables on the 
model response, we examine the behavior of the corresponding first-order 
summands comprising univariate polynomials only, i.e. 

Mf CE { Xi ) = E [M pce (X \Xi = Xi )] (31) 

or equivalently 


= V y, 

aeA 


Ai = {a. e A: OLi> 0, oii^j = 0}. (32) 


Figures [13] and 14 depict such univariate effects considering PCE A and B, 
respectively, for the porosities of the five layers classified as important. De¬ 
spite the different scales, the shapes of respective curves in the two figures 
demonstrate similar trends. The two PCE include univariate polynomials of 
the same degree in <fi C3ab , cj) Lla and <p ci (up to second or third degree), but 
PCE B includes higher-degree univariate polynomials in 0' D4 and (f) Llb (up to 
sixth degree) than PCE A (up to fourth degree). The figures demonstrate 
that an increasing porosity and thus, an increasing hydraulic conductivity 
are associated with a decreasing algebraic contribution to the MLE value. 
Indeed, the hydraulic conductivity parameter, mainly responsible for the ad- 
vective processes within the layer, is linearly and oppositely related to the 
MLE (see Eq. 0 and Q). Overall, the response is more sensitive to changes 






0.8 


Total Sobol’ Indices 




Figure 12: Sobol’ indices using PCE B. 


occurring within low porosity-permeability ranges, where advective processes 
are counterbalanced by dispersive-diffusive processes, with the latter being 
responsible for higher MLE values. Plateaus observed in certain graphs in¬ 
dicate the ascendancy of one ageing process over the other, which results in 
more regular trends. 

For a more in-depth investigation of the contributions of the different 
types of hydro-dispersive parameters, Table [5] lists the sums of the first- 
order indices per type of property over all layers or zones of the considered 
cross-section. According to this table, the added main effects of the porosity 
parameters account for approximately 87% and 93% of the response variance 
for PCE A and B, respectively, whereas the added main effects of all remain¬ 
ing input random variables account for a mere 2.5% and 2.8%, respectively. 
We note the higher contributions of main effects for PCE B and the zero 
main effects of A a for both PCE. 

The above analysis indicated that among the set of random variables, 


29 











Figure 13: Univariate effects of important variables using PCE A. 





Table 5: Sums of first-order Sobol’ indices over all layers per type of property. 


meta-model 

0 

A k 

9 

OtL 

A a 

X7H 

A 

0.8664 

0.0088 

0.0029 

0.0076 

0 

0.0057 

B 

0.9302 

0.0096 

0.0032 

0.0088 

0 

0.0061 


only the porosity parameters of certain layers are important for explaining the 


30 
























response variance. As highlighted earlier, because of the assumed relationship 
between porosity and hydraulic conductivity in each layer, the Sobol’ indices 
for the (j) variables are also indicative of the importance of K x in the respective 
layers. Accordingly, in the subsequent discussion of the SA results, we will 
interpret the variability of MLE in terms of joint effects of the petrofacies P. 


4.4 Discussion of results 


The petrofacies of aquifer formations have the largest effects on the variability 
of the MLE at the TZ. For layers D4 and Lib in particular, for which the 
upper bounds of the hydraulic conductivity ranges are the highest (see Table 
[3} Figure [5}, strong advective processes within their own volume may reduce 
the overall groundwater residence time in the model. Besides, the wider 
the ranges of K x values in these permeable formations, the higher are the 
contributions of the respective petrofacies to the variability of the MLE. The 
three most significant aquifer formations, namely D4, Lib and Lla, have 
rather non-linear univariate effects on the output response (see Figures 13 
and 14). Substantial changes within the response are observed with small 


shifts in the ranges of low porosity-permeability values, revealing the effect of 
the balance between advective and dispersive-diffusive transport processes. 
But for higher porosity-permeability values, advective fluxes prevail, thus 
yielding more linear and moderate changes in the response. 

The position of the aquifers relatively to layer C2 is also relevant: the 
further the layer is, the lower is its effect on the MLE of water molecules 
departing from the TZ. Layer Lla, which is the first aquifer encountered 
in the Oxfordian sequence (at a distance of 60 meters from layer C2), has 
a significant effect on the output variance, whereas layers D2 and D3 that 
have similar K x ranges have much smaller contributions. 

The relatively large uncertainty range and high upper-bound value of per¬ 
meability for layer D1 results in a marginal effect on the response variability 


see Figures 10 and 12), which is nevertheless unimportant with respect to 


the threshold Sj < 0.01. Its remote location relative to layer C2 explains 
the reduced quantity of water molecules departing from the TZ that reach 
the layer, given that most of the solute is captured by the highly advective 
aquifer from Upper Bathonian (D4) situated in between. 

The petrofacies of semi-permeable formations, P C3ab and P cx , are also 
significantly influencing the variability of the MLE at the TZ. By isolating 
layer C2 from major aquifer formations, they buffer solute intrusions into the 
Oxfordian and Dogger aquifer sequences, thus acting like a geological barrier. 


Figures [13] and 14 show that their univariate effects on the output quantity 
are relatively non-linear despite their limited amplitude. 
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We underline that the petrofacies of the host layer attributed to the 
Callovo-Oxfordian claystone ( P c2 ) are insignificant with regard to the vari¬ 
ability of the ouput response. This feature was already outlined in the SA 
performed upon the integrated high-definition model in (Deman et al., 2015). 
Slow diffusive processes take place into highly impermeable rocks, which in¬ 
duces large values of the MLE response (Figure |3| and prove the high effi¬ 
ciency of the Callovo-Oxfordian claystone to act as a geological barrier. Mod¬ 
ifying the magnitudes of advective-dispersive transport processes in layer C2 
does not add a significant variability to the time required for solutes to leave 
the domain where the numerical model is defined. 

The magnitude of the transverse advective fluxes in each layer is related to 
the respective value of K z , the uncertainty of which is accounted for through 
the anisotropy ratio Ak- Although Ak represents the second most influential 
group of factors considering added effects from all layers (see Table [5]), the 
uncertainty in this property adds a small amount of variability to the MLE 
(< 1%) compared to the petrofacies (~ 90%). We note however that fac¬ 
tor A^ 4 is only marginally excluded from the important factors when PCE 
B is considered (see Figure 12). Indeed, for the highest K x values, strong 
advective fluxes take place within the layer’s volume. Under the assump¬ 
tion of strong transverse fluxes (A£- 4 —» 1), solutes can be oriented toward 
neighbouring layers where slower fluxes occur, thus raising the MLE. 

For each layer, the Euler angle 6 could deviate groundwater fluxes from an 
orientation parallel to the x-axis and toward the main discharge boundaries, 
thus raising the variability of the response. Although it could be assumed 
as especially influential in the most advective layers, the total contribution 
of this group of uncertain factors to the variance of the MLE is negligible in 
comparison to that of the petrofacies P (see Table [5]) . 

In aquifer formations, the effects of the uncertainty regarding the macro¬ 
dispersion tensors upon the response quantity, i.e. the magnitude (a^) and 
anisotropy (A a ) of the tensors, are concealed by the strong effect of petrofa¬ 
cies on the advective part of the transport processes (see Table [5]). We note 
however that the longitudinal component of the macro-dispersion tensor in 
layer C3ab (a£ 3ab ) appears among the ten factors with the highest total 
Sobol’ indices for both PCE A and B. The anisotropy ratios, A a , have no 
contribution at all to the response variance when considered independently; 
the uncertainty in these factors contributes to the variability of MLE only 
through interaction terms. 

The sensitivity of the MLE with respect to flow BCs considered in the 
model is directly related to the magnitude and orientation of the advective 
fluxes in the entire model. In the case of high gradients in both limestone 
sequences, the advective solute transport processes would raise within their 
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volume, and thus reduce the MLE. Table [5] indicates that the three random 
hydraulic gradients have a small added effect on the output variance. Note 
however that the total Sobol’ index for the hydraulic gradient in the Dogger 
sequence (V// 1 ) belongs to the ten highest indices for both PCE A and 
B. The uncertainty regarding this factor can alter the advective processes 
occurring notably in layer D4, which has the far highest contribution to the 
variance of the MLE calculated at the TZ. 


5 Conclusions 

The model introduced in this paper stands as a two-dimensional vertical 
cross-section of the subsurface of Paris Basin in the vicinity of Bure (Haute- 
Marne). While encompassing most of the hydrogeological features of the 
underground media, it was simplified with regard to geometries, disconti¬ 
nuities, fractures and heterogeneities. This numerical model is intended to 
explore the behavior of a complex multi-layered hydrogeological system at 
low computational cost and provide insights into the effect of uncertain pa¬ 
rameters upon the mean lifetime expectancy (MLE). 

Sensitivity analysis (SA) was carried out considering a high-dimensional 
random input. For the sake of simplicity, homogeneous parameters were 
assumed within each of the 15 hydrogeological layers comprising the model. 
The uncertain factors at each layer included: the petrofacies, P, regarded 
as the couple permeability-porosity, {K, (j)}; the anisotropy ratio and the 
orientation of the hydraulic conductivity tensor, Ax and 6 respectively; the 
magnitude and anisotropy ratio in the macro-dispersion tensor, (Xl and A a 
respectively. Additionally, the hydraulic gradients, VH, in three zones of the 
model domain were considered random, leading to a total of 78 uncertain 
input factors. 

In the present study, a target zone (TZ) located within the middle layer 
(C2) of the model domain provides the output response of interest. Latin 
hypercube sampling was employed to address the propagation of the uncer¬ 
tainty from the input factors upon the MLE of water molecules departing 
from the TZ. Polynomial chaos expansion (PCE) meta-models were used to 
compute the Sobol’ sensitivity indices for each input factor at low compu¬ 
tational costs. Sparse PCE proved singularly efficient in providing accurate 
representations of the response of interest at low computational cost with 
respect to the high dimensionality of the model. The accuracy was enhanced 
when the PCE were fitted to the logarithmic MLE; because of the wide range 
of variation of the MLE, considering its logarithmic transformation is herein 
meaningful from a physical standpoint. 
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The SA results for the non-transformed and the logarithmic MLE demon¬ 
strated similar trends. It was found that the variability of the MLE is almost 
entirely due to the uncertainty regarding the petrofacies of the hydrogeolog¬ 
ical layers. The other hydro-dispersive parameters are insignificant for ex¬ 
plaining the response variance and may be considered as deterministic factors 
in future works. 

Focusing on the effects of petrofacies solely, the SA demonstrated the large 
contributions of aquifer formations to the variance of the model output. In 
particular, (i) the closer the aquifer formation to layer C2, (ii) the thicker 
the layer, (iii) the wider the ranges of the petrofacies and (iv) the higher 
the upper bound of the hydraulic conductivity range, the larger were the 
effects of the petrofacies on the variability of the response. Investigation 
of the univariate effects of petrofacies highlighted that for these permeable 
formations, the response is more sensitive to changes occurring within low 
porosity-permeability ranges. Hence, within a certain range of {A', 0} values, 
the dispersive-diffusive processes counterbalance with the strong advective 
fluxes in the ageing process. 

In formations characterized by highly advective processes, the longitudi¬ 
nal hydraulic conductivities applying in the main groundwater direction have 
large contributions to the MLE variability. The two semi-confining forma¬ 
tions encompassing the C2 layer buffer the transfer of solute from the latter 
toward the further aquifer sequences. Besides, it is acknowledged that lon¬ 
gitudinal dispersion processes occurring within their own volume also retard 
the solute transfer toward the adjacent aquifers. Because of the diffusion- 
dominated transport processes occurring within its volume, the petrofacies 
of the highly-confining C2 layer have a negligible effect on the variance of 
the output response although responsible for the high values of the latter. 

It is important to remind that the use of a 2D model tends to underesti¬ 
mate the output response of interest by omitting the advection and dispersion 
along the third dimension. Recognizing this limitation, we underline that the 
purpose of the simplified model introduced herein is to shed light on the rela¬ 
tive effects of various uncertain factors governing the advective and dispersive 
processes in a complex multi-layered hydrogeological system. The presented 
methodology may be applied to a real-case study employing a realistic 3D 
numerical model. 

The sensitivity analysis performed in this work is deemed particularly in¬ 
formative for future applications with the high-resolution integrated Meuse/Haute- 
Marne hydrogeological model. In the frame of a real-case uncertainty analysis 
with concern to a solute transport in the subsurface of the Paris Basin, the 
authors recommend defining as thoroughly as possible the spatial distribu¬ 
tions of hydraulic conductivity values, with a main focus on the large aquifer 
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sequences of Oxfordian and Dogger ages closer to the potential host layer. 


Acknowledgement 

The authors are grateful to Mr Benjamin Brigand (Universite Paris-Sud) 
and Mrs Agnes Vinsot (ANDRA) for providing the hydro-dispersive datasets 
used in this analysis, and to Mr Fabien Cornaton (DHI-WASY GmbH) for 
providing the code GroundWater. 


References 

Renard, P.. Stochastic hydrogeology: What professionals really need? 
Ground Water 2007;45(5):531 541. 

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, 
D., et al. Global Sensitivity Analysis - The Primer. Wiley; 2008. 

Cukier, H., Levine, R., Shuler, K.. Nonlinear sensitivity analysis of multi¬ 
parameter model systems. J Comput Phys 1978;26:1-42. 

Saltelli, A., Tarantola, S., Chan, K.. A quantitative, model independent 
method for global sensitivity analysis of model output. Technometrics 
1999;41 (1) :39 56. 

Sobol’, I.M.. Sensitivity estimates for nonlinear mathematical models. Math 
Modeling & Comp Exp 1993;1:407-414. 

Archer, G., Saltelli, A., Sobol’, I.. Sensitivity measures, ANOVA-like 
techniques and the use of bootstrap. J Stat Comput Simul 1997;58:99- 
120 . 

Sobol’, I.. Global sensitivity indices for nonlinear mathematical models and 
their Monte Carlo estimates. Mathematics and Computers in Simulation 
2001 ;55 (1-3) :271-280. 

Saltelli, A.. Making best use of model evaluations to compute sensitivity 
indices. Comput Phys Comm 2002;145:280-297. 

Sobol’, I., Kucherenko, S.. Global sensitivity indices for nonlinear mathe¬ 
matical models. Review. Wilmott magazine 2005;1:56-61. 

Saltelli, A., Annoni, P., Azzini, V., Campolongo, F., Ratto, M., Tarantola, 
S.. Variance based sensitivity analysis of model output. Design and estima¬ 
tor for the total sensitivity index. Comput Phys Comm 2010;181:259-270. 


35 



Sobol, I.M., Tarantola, S., Gatelli, D., Kucherenko, S.S., Mauntz, 
W.. Estimating the approximation error when hxing unessential factors 
in global sensitivity analysis. Reliability Engineering & System Safety 
2007;92(7):957-960. 

Janon, A., Klein, T., Lagnoux, A., Nodet, M., Prienr, C.. Asymptotic nor¬ 
mality and efficiency of two Sobol index estimators. ESAIM: Probability 
and Statistics 2013;: 1 20. 

Oakley, J., O’Hagan, A.. Probabilistic sensitivity analysis of complex 
models: a Bayesian approach. J Royal Stat Soc, Series B 2004;66:751-769. 

Marrel, A., Iooss, B., Laurent, B., Roustant, O.. Calculations 

of Sobol indices for the Gaussian process metamodel. Rcliab Eng 
Sys Safety 2009;94:742-751. URL http://hal.archives-ouvertes.fr/ 
hal-00239494/PDF/SAM0_marrel_hal.pdf, 

Storlie, C.B., Swilcr, L.P., Helton, J.C., Sallaberry, C.J.. Implementa¬ 
tion and evaluation of nonparametric regression procedures for sensitivity 
analysis of computationally demanding models. Reliab Eng Sys Safety 
2009;94:1735-1763. 

Zuniga, M.M., Kucherenko, S., Shah, N.. Metamodelling with in¬ 
dependent and dependent inputs. Computer Physics Communications 
2013; 184(6): 1570-1580. 

Sudret, B.. Global sensitivity analysis using polynomial chaos expansions. 
Rehab Eng Sys Safety 2008;93:964-979. 

Blatman, G., Sudret, B.. Efficient computation of global sensitivity in¬ 
dices using sparse polynomial chaos expansions. Reliab Eng Sys Safety 
2010a;95:1216-1229. 

Fajraoui, N., Ramasomanana, F., Younes, A., Mara, T., Ackerer, 
P., Guadagnini, A.. Use of global sensitivity analysis and polynomial 
chaos expansion for interpretation of nonreactive transport experiments in 
laboratory-scale porous media. Water Resources Research 2011;47. 

Younes, A., Mara, T., Fajraoui, N., Lehmann, F., Belfort, B., Bey- 
doun, H.. Use of global sensitivity analysis to help assess unsaturated soil 
hydraulic parameters. Vadose Zone Journal 2013; 12(1). 

Sochala, P., Le Maitre, O.. Polynomial chaos expansion for subsurface flows 
with uncertain soil parameters. Advances in Water Resources 2013;62:139- 
154. 


36 



Ciriello, V., Federico, V., Riva, M., Cadini, F., Sanctis, J., Zio, E., 
et al. Polynomial chaos expansion for global sensitivity analysis applied 
to a model of radionuclide migration in a randomly heterogeneous aquifer. 
Stochastic Environmental Research and Risk Assessment 2013;27(4) :945 
954. 

Formaggia, L., Guadagnini, A., Imperiali, I., Lever, V., Porta, G., Riva, 
M., et al. Global sensitivity analysis through polynomial chaos expansion of 
a basin-scale geochemical compaction model. Computational Geosciences 
2013;17(l):25-42. 

Laloy, E., Rogiers, B., Vrugt, J.A., Mallants, D., Jacques, D.. Effi¬ 
cient posterior exploration of a high-dimensional groundwater model from 
two-stage markov chain monte carlo simulation and polynomial chaos ex¬ 
pansion. Water Resources Research 2013;49(5):2664-2682. 

Delay, J., Trouillcr, A., Lavanchy, J.. Hydrodynamic properties of the 
Callovo-Oxfordian formation in the East of the Paris Basin: comparison of 
results obtained through different approaches. Comptes Rendus Geoscience 
2006;338(12-13) :892 907. 

Distinguin, M., Lavanchy, J.. Determination of hydraulic properties of 
the Callovo-Oxfordian argillite at the Bure site: Synthesis of the results 
obtained in deep boreholes using several in situ investigation techniques. 
Physics and Chemistry of the Earth 2007;32(l-7):379-392. 

Enssle, C., Cruchaudet, M., Croise, J., Brommundt, J.. Determination of 
the permeability of the Callovo-Oxfordian clay at the metre to decametre 
scale. Physics and Chemistry of the Earth 2011;36(17-18) :1669 1678. 

Brigand, B., Vincent, B., Durlct, C., Deconinck, J., Blanc, P., Trouiller, 
A.. Acoustic properties of ancient shallow-marine carbonates: Effects of 
depositional environments and diagenetic processes (middle jurassic, Paris 
Basin, France). Journal of Sedimentary Research 2010;80(9-10):791-807. 

Linard, Y., Vinsot, A., Vincent, B., Delay, J., Wechner, S., 
De La Vaissiere, R., et al. Water flow in the Oxfordian and Dogger lime¬ 
stone around the Meuse/Haute-Marne underground research laboratory. 
Physics and Chemistry of the Earth 2011;36(17-18):1450-1468. 

Landrein, P., Vigneron, G., Delay, J., Lebon, P., Pagel, M.. Lithology, hy¬ 
drodynamism and thermicity in the multi-layer sedimentary system inter¬ 
sected by the Andra deep borehole of Montiers-sur-Saulx (Meuse, France). 
Bulletin de la Societe Geologique de France 2013;184(6):519-543. 


37 



Deman, G., Kerrou, J., Benabderrahmane, H., Perrochet, P.. Sensitivity 
analysis of groundwater lifetime expectancy to hydro-dispersive parame¬ 
ters: the case of ANDRA Meuse/Haute-Marne site. Reliability Engineering 
& System Safety 2015; 134(0):276-286. 

ANDRA. Modelc hydrogeologique integre region-secteur a l’actuel. Phase 
11 2010-2012 - tache 3 : Sensibilite du modele de reference aux incer¬ 
titudes sur les parametres hydro-dispersifs et analyse de risque. Report 
C.RP.0CHYN.12.0001. Tech. Rep.; 2012. 

Castellini, A., Chawathe, A., Larue, D., Landa, J., Jian, F., Toldi, J., 
et al. What is relevant to flow? a comprehensive study using a shallow 
marine reservoir. 2003. 

Bourgeat, A., Kern, M., Schumacher, S., Talandier, J.. The COUPLEX 
test cases: Nuclear waste disposal simulation. Computational Geosciences 
2004;8(2):83-98. 

Uffink, G.. Application of the Kolmogorov’s backward equation in ran¬ 
dom walk simulation of groundwater contaminant transport. Contaminant 
Transport in Groundwater 1989; :283 289. 

Cornaton, F., Perrochet, P.. Groundwater age, life expectancy and transit 
time distributions in advective-dispersive systems; 2. reservoir theory for 
sub-drainage basins. Advances in Water Resources 2006a;29(9):1292-1305. 

Cornaton, F., Perrochet, P.. Groundwater age, life expectancy and transit 
time distributions in advective-dispersive systems: 1. generalized reservoir 
theory. Advances in Water Resources 2006b;29(9):1267-1291. 

Kazemi, G., Lehr, J., Perrochet, P.. Groundwater Age. Wiley; 2006. 

Cornaton, F.. GroundWater: a 3-D ground water and surface water flow, 
mass transport and heat transfer finite element simulator, reference man¬ 
ual. 2007. 

Kerrou, J., Renard, P.. A numerical analysis of dimensionality and hetero¬ 
geneity effects on advective dispersive seawater intrusion processes. Hy¬ 
drogeology Journal 2010;18(1) :55 72. 

Fourre, E., Jean-Baptiste, P., Dapoigny, A., Laviclle, B., Smith, T., 
Thomas, B., et al. Dissolved helium distribution in the Oxfordian and 
Dogger deep aquifers of the Meuse/Haute-Marne area. Physics and Chem¬ 
istry of the Earth 2011;36(17-18):1511—1520. 


38 



Delay, J., Distinguin, M.. Hydrogeological investigations in deep wells 
at the Meuse/Haute Marne underground research laboratory; vol. 104 of 
Lecture Notes in Earth Sciences ; chap. 26. Springer Berlin Heidelberg; 
2004, p. 219-225. 

Delay, J., Distinguin, M., Dewonck, S.. Characterization of a clay-rich 
rock through development and installation of specific hydrogeological and 
diffusion test equipment in deep boreholes. Physics and Chemistry of the 
Earth 2007a;32(l-7):393-407. 

Cosenza, P., Ghoreychi, M., de Marsily, G., Vasseur, G., Violette, S.. 
Theoretical prediction of poroelastic properties of argillaceous rocks from 
in situ specific storage coefficient. Water Resources Research 2002;38(10). 

Delay, J., Rebours, H., Vinsot, A., Robin, P.. Scientific investigation in 
deep wells for nuclear waste disposal studies at the Meuse/Haute Marne 
underground research laboratory, Northeastern France. Physics and Chem¬ 
istry of the Earth 2007b;32(l-7):42-57. 

Mazurek, M., Alt-Epping, P., Bath, A., Ginnni, T., W.H., N., Buschaert, 
S., et al. Natural tracer profiles across argillaceous formations. Applied 
Geochemistry 2011;26(7): 1035-1064. 

Vinsot, A., Delay, J., de La Vaissiere, R., Cruchaudet, M.. Pumping tests 
in a low permeability rock: Results and interpretation of a four-year long 
monitoring of water production flow rates in the Callovo-Oxfordian argilla¬ 
ceous rock. Physics and Chemistry of the Earth, Parts A/B/C 2011;36(17- 
18): 1679-1687. 

Contoux, C., Violette, S., Vivona, R., Goblet, P., Patriarche, D.. How 
basin model results enable the study of multi-layer aquifer response to 
pumping: the Paris Basin, France. Hydrogeology Journal 2013;21(3):545- 
557. 

de Hoyos, A., Viennot, P., Ledoux, E., Matray, J., Rocher, M., Certes, C.. 
Influence of thermohaline effects on groundwater modelling - application 
to the Paris sedimentary Basin. Journal of Hydrology 2012;464:12-26. 

Goncalves, J., Violette, S., Guillocheau, F., Robin, C., Pagel, M., Bruch 
D., et ah Contribution of a three-dimensional regional scale basin model 
to the study of the past fluid flow evolution and the present hydrology of 
the Paris Basin, France. Basin Research 2004a;16(4):569-586. 


39 



Goncalves, J., Violette, S., Robin, C., Bruel, D., Guillocheau, F., Ledoux, 
E.. Combining a compaction model with a facies model to reproduce 
permeability fields at the regional scale. Physics and Chemistry of the 
Earth 2004b;29(l):17-24. 

Li, G., Rabitz, H., Yelvington, P., Oluwole, O., Bacon, F., Kolb, C., et al. 
Global sensitivity analysis for systems with independent and/or correlated 
inputs. J PhysChem 2010;114:6022-6032. 

Kucherenko, S., Tarantola, A., Armorii, P.. Estimation of global sensi¬ 
tivity indices for models with dependent variables. Comput Phys Comm 
2012;183:937-946. 

Mara, T.A., Tarantola, S.. Variance-based sensitivity indices for mod¬ 
els with dependent inputs. Reliability Engineering & System Safety 
2012;107:115-121. 

Xiu, D., Karniadakis, G.E.. The Wiener-Askey polynomial chaos for 
stochastic differential equations. SIAM J Sci Comput 2002;24(2):619-644. 

Blatman, G., Sudret, B.. An adaptive algorithm to build up sparse poly¬ 
nomial chaos expansions for stochastic finite element analysis. Prob Eng 
Mech 2010b;25:183-197. 

Tibshirani, R... Regression shrinkage and selection via the Lasso. J Royal 
Stat Soc, Series B 1996;58:267-288. 

Blatman, G., Sudret, B.. Adaptive sparse polynomial chaos expansion 
based on Least Angle Regression. J Comput Phys 2011;230:2345-2367. 

Efron, B., Hastie, T., Johnstone, I., Tibshirani, R.. Least angle regression. 
Annals of Statistics 2004;32:407-499. 

Allen, D.. The prediction sum of squares as a criterion for selecting predictor 
variables. Technical report, Dept of Statistics, LIniversity of Kentucky 
1971; (23). 

Chapellc, O., Vapnik, V., Bengio, Y.. Model selection for small sample 
regression. Machine Learning 2002;48(l):9-23. 

Marelli, S., Sudret, B.. UQLab: a framework for uncertainty quantification 
in MATLAB. In: Proc. 2nd Int. Conf. on Vulnerability, Risk Analysis and 
Management (ICVRAM2014), Liverpool, LInited Kingdom. 2014,. 


40 



Marclli, S., Lamas-Fernandes, C., Sudret, B.. Uqlab user manual - sensi¬ 
tivity analysis. Tech. Rep.; Chair of Risk, Safety & Uncertainty Quantifi¬ 
cation, ETH Zurich; 2015. Report # UQLab-VO.9-106. 

McKay, M., Beckman, R., Conover, W.. A comparison of three methods 
for selecting values of input variables in the analysis of output from a 
computer code. Technometrics 19T9;21 (2) :239 245. 

Wang, G.. Adaptive response surface method using inherited Latin Hyper¬ 
cube design points. J Mech Design 2003;125:210-220. 

Blatman, G., Sudret, B.. An adaptive algorithm to build up sparse poly¬ 
nomial chaos expansions for stochastic finite element analysis. Prob Eng 
Mech 2010c;25:183-197. 

Borgonovo, E., Tarantola, S., Plischke, E., Morris, M.. Transformations 
and invariance in the sensitivity analysis of computer experiments. Jour¬ 
nal of the Royal Statistical Society: Series B (Statistical Methodology) 
2014;76(5):925-947. 


41 



